Problème torsion d’arbres

Problème torsion d’arbres#

Librairies#

# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')

Géométrie#

Définition de la géométrie et du maillage#

import gmsh
import numpy as np
import meshio
from dolfinx.io import XDMFFile
from mpi4py import MPI
import ipywidgets as widgets
from IPython.display import display
import pyvista as pv
from dolfinx import plot

# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")

# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
    global domain
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("ellipse")

    # Création de l'ellipse
    ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
    curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
    surface = gmsh.model.occ.addPlaneSurface([curve_loop])

    # Extrusion pour obtenir un cylindre elliptique
    gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
    gmsh.model.occ.synchronize()

    # Configuration du maillage
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
    gmsh.model.mesh.generate(3)

    # Sauvegarde et conversion du maillage
    gmsh.write("ellipse.msh")
    msh = meshio.read("ellipse.msh")
    meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
    gmsh.finalize()

    # Lecture du maillage avec FEniCSx
    with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
        domain = xdmf_file.read_mesh(name="Grid")
        domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
    
    # Conversion du maillage pour la visualisation avec PyVista
    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
    u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

    # Effacer la scène précédente et ajouter le nouveau maillage
    p.clear()
    # Ajout du maillage à la scène de visualisation
    p.add_mesh(u_grid, 
               show_edges=True,  # Affiche les arêtes du maillage
               scalar_bar_args={
                   "title": "u",  # Titre de la barre de couleur
                   "title_font_size": 24,
                   "label_font_size": 22,
                   "shadow": True,
                   "italic": True,
                   "font_family": "arial",
                   "vertical": False  # Orientation horizontale de la barre de couleur
               })

    # Ajout d'un titre à la visualisation
    p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")

    # Ajout des limites de la boîte englobante
    p.show_bounds(color="black")

    # Ajout des axes de coordonnées
    p.add_axes(color="black")

    # Définition de la couleur de fond
    p.set_background("grey")
    
    # Mise à jour de la scène
    p.show()
    
    return rho_x, rho_y, h
    
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')

# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})

display(ui, out)

# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
    current_rho_x = rho_x_slider.value
    current_rho_y = rho_y_slider.value
    current_h = h_slider.value
    return current_rho_x, current_rho_y, current_h

# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y])  # Demi-axes dans les directions x et y

Configuration du problème#

Définition de l’espace fonctionnel#

# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Définition des frontière du domaine#

Visualisation des frontière du domaine#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites#

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau#

import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox

# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")

lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")

# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])

# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
    mu = mu_input.value * 1e9  # Conversion en Pa
    lambda_ = lambda_input.value * 1e9  # Conversion en Pa
    
    # Afficher les valeurs en Pa
    clear_output(wait=True)  # Clear previous output
    
    return mu, lambda_  # Retourner les valeurs en Pa

# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')

# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)

mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")

Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00

Définition des chargements#

import ipywidgets as widgets
from IPython.display import display
import numpy as np

def create_angle_widget():
    # Création du widget curseur
    angle_slider = widgets.FloatSlider(
        value=20.0,  # Valeur initiale à 20 degrés
        min=0,
        max=360,
        step=1,
        description='Angle (°):',
        disabled=False,
        continuous_update=True,
        orientation='horizontal',
        readout=True,
        readout_format='.0f',
    )

    # Widget de sortie pour afficher les valeurs
    output = widgets.Output()

    # Fonction pour mettre à jour l'affichage
    def update_angle(change):
        ang_degrees = change['new']
        ang_radians = np.radians(ang_degrees)
        
        with output:
            output.clear_output(wait=True)
            print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
        
        return ang_radians

    # Attacher la fonction de mise à jour au widget
    angle_slider.observe(update_angle, names='value')

    # Afficher le widget et la sortie
    display(widgets.VBox([angle_slider, output]))

    # Initialiser l'affichage
    update_angle({'new': angle_slider.value})

    return angle_slider

# Créer et afficher le widget
angle_widget = create_angle_widget()

# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
    return np.radians(angle_widget.value)

Taux de torsion (α) = 3.0020 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 6.0358e+08 N·m²

# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4)  # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")

# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)

# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1]  # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0]   # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0)  # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)

# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")

Densité de contrainte de torsion (A_τ) = 2.4016e+11 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)

Résolution#

# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
    options=[
        ('Aucune', 'none'),
        ('Face inférieure S0', 'S0'),
        ('Face supérieure Sh', 'Sh'),
        ('Surface latérale Sl', 'Sl')
    ],
    value='S0',  # Valeur par défaut
    description='Encastrement:',
    disabled=False,
)

# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
    'none': [],
    'S0': [bc0_S0],
    'Sh': [bc0_Sh],
    'Sl': [bc0_Sl]
}

# Fonction pour mettre à jour le problème
def update_problem(change):
    selected_bc = change['new']
    problem = LinearProblem(a, L, bcs=bc_dict[selected_bc], 
                            petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
    
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Encastrement appliquée sur {selected_bc}")
            
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')

# Afficher le widget
display(bc_selector, output)

# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})

uh = problem.solve()

Visualisation des résultats#

Visualisation des déplacements#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats#

Calcul des valeurs maximales#

Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 3.001966 radians
Angle de torsion théorique θ : 172.00 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 2.985258 radians
Rotation Z maximale simulée : 171.04 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 2.987141 radians
Norme de rotation maximale : 171.15 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 18.72 degrés
RX max  : 0.326659 radians
RY max : 18.33 degrés
RY max  : 0.320002 radians

Problème torsion d’arbres

Contents

Problème torsion d’arbres#

Librairies#

Géométrie#

Définition de la géométrie et du maillage#

Configuration du problème#

Définition de l’espace fonctionnel#

Définition des frontière du domaine#

Visualisation des frontière du domaine#

Conditions aux limites#

Propriétés matériau#

Définition des chargements#

Résolution#

Visualisation des résultats#

Visualisation des déplacements#

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

Visualisation des contraintes de von Mises#

Analyses des résultats#

Calcul des valeurs maximales#